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Introduction 

This  document  reports  the  progress  in  the  third  year  of  a  four  year  grant.  This  introduction 
includes  a  brief  description  of  the  subject,  goals,  and  scope  of  the  research  project,  along  with  a 
description  of  the  relationship  of  this  work  with  that  of  previous  researchers.  Following  this 
description  is  a  milestone  chart  and  an  outline  of  the  progress  to  date.  (This  would  be  a  good  place 
to  start  for  a  brief  overview  of  the  current  program  status.)  The  introductory  section  ends  with  a 
more  detailed  review  of  the  results  of  the  first  two  years  of  the  program.  The  main  body  of  the 
report  is  a  detailed  description  of  the  accomplishments  during  the  third  year  of  the  program.  This 
is  followed  by  some  concluding  remarks  summing  up  the  current  state  of  the  research  and  future 
tasks  to  be  accomplished. 

Program  Description  and  Goals 

Hyperthermia  has  been  shown  to  be  an  efficacious  adjuvant  therapy  in  the  treatment  of 
recurrent  breast  cancer.  [1,2,3]  One  of  the  most  flexible  and  proven  methods  to  induce 
hyperthermia  is  ultrasound.  In  order  to  optimally  utilize  such  clinical  devices  it  is  essential  that  the 
operator  have  the  ability  to  accurately  plan  a  treatment.  The  clinical  hyperthermia  situation  is 
diagrammed  in  Figure  1.  A  treatment  plan  consists  of  a  defined  motion  path  and  applied  transducer 
power  that  will  raise  the  temperature  of  the  tumor  to  a  desired  set  point  for  a  specified  time  while 
minimizing  the  effect  on  surrounding  tissue.  The  treatment  plan  must  be  adjusted  for  tissue 
geometry,  absorption,  and  desired  temperature  history  specific  to  each  patient.  Planning  a 
treatment  begins  with  the  accurate 
estimation  of  the  absorbed  power  field 
created  by  the  ultrasound  transducer. 

Thus,  the  ability  to  accurately  estimate 
this  field  is  critical  to  producing  a  proper 
and  effective  treatment  plan. 

It  is  the  goal  of  this  research  to 
provide  improved  numerical  modeling  of 
ultrasound  propagation  through  breast 
tissue  and  the  post-mastectomy  chest  wall 
for  use  in  patient  treatment  planning  of 
hyperthermia  cancer  treatments.  With  this 
improved  model  the  clinician  can  plan  an 
ultrasound  hyperthermia  treatment  so  as 
to  produce  the  required  thermal  dose 
while  minimizing  patient  pain  due  to 
excessive  temperature  or  ultrasound- 
tissue  interactions.  A  clinical  system  is 
being  constructed  to  combine  model 
results  with  experiments  in  an  inverse 
technique  to  estimate  the  absorbed  power 
field  in  the  breast  and  chest  wall  during 
ultrasound  hyperthermia.  This  improved 
planning  will  aid  the  clinicians  in 
providing  better  hyperthermia  treatments, 
which  will  improve  treatment  response 
rates. 

The  system  under  development  consists  of  the  following  subsystems: 

1 .  An  ultrasonic  B-mode  imager  based  clinical  data  acquisition  system  used  to  obtain  patient 

anatomy  and  measure  the  attenuation  of  ultrasound  by  the  body  tissue. 


Patient 


Tumor 


'Membrane 


Water  Bath 


Ultrasound 

Propagation 


Ultrasound  Transducer 

I  Figure  1 :  Diagram  of  the  clinical  situation  for  a  hyperthermia 
treatment  of  the  breast  or  chest  wall  area. 
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2.  A  clinical  technique  using  an  ultrasound  treatment  system  to  directly  measure  absorbed 
power  at  specific  points. 

3.  An  improved  model  of  ultrasound  propagation  through  breast  tissue. 

4.  An  inverse  technique  that  integrates  experimental  measurements  and  modeling  results  to 
obtain  the  “best"  prediction  of  the  ultrasound  power  deposition  field. 

The  B-mode  diagnostic  ultrasound  imager  in  Subsystem  1  is  used  to  construct  a  geometric 
model  based  upon  each  patient's  anatomy  and  to  measure  the  ultrasound  attenuation  coefficient  at 
selected  areas  within  the  tissue  region.  In  Subsystem  2  the  clinical  treatment  system  is  used  in  a 
scanning  protocol  to  locate  thermocouples  imbedded  in  the  breast  tissue.  The  response  of  these 
thermocouples  to  the  scanned  ultrasound  gives  a  direct  measure  of  the  power  deposition  at  that 
point.  Subsystem  3  is  an  ultrasound  propagation  and  power  deposition  model.  For  this  model  we 
have  considered  2  options: 

a)  Hybrid  Green’ s  Function  Model  -  Ultrasound  propagation  from  the  transducer  into  the 
breast  is  modeled  with  a  hybrid  of  Green's  function  solutions  used  in  the  homogenous 
water  region  and  a  finite  element  solution  of  the  acoustic  wave  equation  in  the 
inhomogeneous  breast  region. 

b)  Parabolic  Model  -  This  is  a  forward  marching  parabolic  model  that  is  quite  fast,  but  does 
not  take  into  account  reflected  energy  at  interfaces. 

Subsystem  4  is  an  iterative  inverse  technique  that  optimizes  the  parameters  of  the  power  deposition 
model  by  forcing  the  model  results  to  match  the  experimental  measurements.  The  model  is  used  to 
predict  the  power  deposition  at  the  thermocouple  locations,  the  error  between  model  and 
experiment  is  measured,  and  the  error  values  are  used  to  adjust  the  model  parameters  for  the  next 
iteration. 

This  improved  power  deposition  model  will  be  useful  in  the  development  of  more  accurate 
path  planning  for  hyperthermia  treatment.  As  such  it  will  provide  a  clinically  useful  tool  that  can 
serve  as  the  base  for  the  development  of  more  accurate  treatment  path  planning,  more  accurate 
control  techniques,  and  eventually,  higher  success  rates  in  hyperthermia  treatment  with  reduced 
patient  discomfort  during  treatment. 


Task  Description  and  Milestones: 
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Outline  of  Accomplishments: 

1st  Year  Accomplishments: 

Subsystem  1 :  Geometry  Acquisition  and  Attenuation  Measurement 

•  Acquisition  of  ultrasonic  B-scan  instrument. 

•  Development  of  interface  hardware  to  isolate  analog  A-line  information. 

•  Acquisition  of  high  speed  A/D  system  and  development  of  control  software. 

•  Investigation  of  reflection  based  acoustic  attenuation  measurement  techniques. 

Subsystem  2:  Experimental  SAR  Measurement 

•  Initial  discussions  and  verification  of  equipment. 

Subsystem  3:  Model  Development 

•  Development  and  testing  of  a  2D  finite  element  based  model. 

•  Identification  of  alternative  solutions. 

Subsystem  4:  Inverse  Technique 

•  Initial  discussions. 

2nd  Year  Accomplishments: 

Subsystem  1 :  Geometry  Acquisition  and  Attenuation  Measurement 

•  Development  of  Log  Spectral  Difference  software  for  measuring  the  acoustic 
attenuation  from  a  single  A-line. 

•  Development  and  testing  of  a  Force/Balance  test  apparatus  for  direct  measurement  of 
acoustic  attenuation. 

Subsystem  2:  Experimental  SAR  Measurement 

•  No  work  accomplished. 

Subsystem  3 :  Model  Development 

•  Extension  of  the  2D  FE  model  to  three  dimensions.  Combination  of  the  FE  model 
with  an  integral  method  model  to  improve  the  speed  of  the  model. 

•  Development  of  a  3D  parabolic  model  as  a  much  faster,  possibly  less  accurate 
alternative. 

•  Development  of  a  3D  Green’s  Function  model  for  use  as  the  “gold  standard”  with 
which  the  other  models  will  be  compared. 

•  Initial  model  comparisons. 

Subsystem  4:  Inverse  Technique 

•  No  progress. 

3rd  Year  Accomplishments: 

Subsystem  1:  Geometry  Acquisition  and  Attenuation  Measurement 

•  Laboratory  hardware  was  developed  to  provide  3D  scanning  capability  for  the  medical 
B-scan  imager.  Clinical  hardware  is  now  essentially  complete. 

•  Software  to  produce  a  3D,  tissue  attenuation  map  from  the  A-mode  signals  was 
completed  and  is  undergoing  verification  testing. 
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•  Software  and  hardware  for  clinical  patient  geometry  acquisition  is  approximately  50% 
complete. 

Subsystem  2:  Experimental  SAR  Measurement 

•  Thermocouple  based  beam  plots  have  been  performed  to  calibrate  the  ultrasound 
power  field  in  water  for  specific  transducers. 

•  Experiments  have  been  conducted  in  water  and  are  beginning  using  cast  agar 
phantoms. 

•  Preliminary  in-vivo  experiments  were  conducted  on  a  dog. 

•  A  finite  difference  model  has  been  constructed  to  help  interpret  the  experimental 
measurements. 

Subsystem  3:  Model  Development 

•  The  parabolic  model  has  been  selected  for  use  in  regions  of  low  acoustic  contrast.  In 
these  regions  the  model  is  quite  fast  and  loses  little  accuracy  when  compared  to  the 
Green’s  Function  Model.  The  Green’s  Function  Model  can  be  used  in  high  contrast 
areas  (close  to  ribs,  etc.) 

•  An  interface  has  been  built  for  the  model  which  allows  the  user  to  specify  the  volume 
modeled,  the  acoustic  power  input  as  a  2D  complex  pressure  function,  and  the  model 
parameters;  density,  sound  velocity,  and  absorption  at  each  node  within  the  volume. 

Subsystem  4:  Inverse  Technique 

•  Basic  algorithm  development  is  complete.  Optimization  runs  have  been  completed  on 
simple  geometries. 


First  Year  Summary 

The  first  year  efforts  centered  on  initial  development  of  the  ultrasound  propagation  models 
and  acquisition  and  initial  development  of  the  diagnostic  ultrasound  equipment.  A  brief  discussion 
of  these  efforts  is  included  below.  For  more  detailed  information  refer  to  the  first  year  progress 
report. 

Initial  modeling  efforts  focused  on  two  potential  methods:  a  Green’s  Function  Method 
where  the  transducer  face  is  modeled  as  a  distribution  of  point  sources,  and  the  Finite  Element 
Method  where  FE  techniques  are  used  to  solve  the  acoustic  wave  equation.  A  two  dimensional 
FE  code  was  written,  and  verification  runs  were  under  way  at  the  end  of  the  first  year.  The 
results  of  these  runs  and  further  investigation  of  the  modeling  possibilities  led  to  the  development 
of  a  hybrid  model  that  combines  the  features  of  the  two  techniques  to  provide  improved  speed 
and  accuracy.  Details  of  this  model  are  discussed  in  the  second  year  accomplishments  below. 

During  the  first  year  a  diagnostic  ultrasound  B-mode  imager  was  acquired.  This  imager 
was  modified  to  allow  access  to  the  raw  RF  waveform  (A-line)  for  each  scan  line  of  the  B-mode 
image.  A  schematic  of  the  resulting  system  is  shown  in  Figure  2.  The  ultrasound  unit  (Scanner 
250,  Pie  Medical  USA  3535  Route  66,  Neptune  NJ  07753)  was  modified  to  output  the  RF 
signals  along  with  sync  signals.  A  custom  hardware  selection  device  was  constructed  to  read 
these  signals  and  output  a  single  A-line.  A  PC  based  data  acquisition  system  was  acquired  that 
allows  digitization  of  the  A-line  signals.  Once  the  data  is  in  the  PC,  it  is  available  for  use  in 
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software  techniques  to  estimate  the  ultrasonic  attenuation  field  in  the  tissue.  Along  with  the 
hardware  development,  software  was  developed  to  operate  the  hardware  and  perform  the 
attenuation  estimate.  [4] 

Second  Year  Accomplishments 

Subsystem  1:  Ultrasound  Based  Anatomy  and  Attenuation  Measurement 

Ultrasonic  Attenuation  Measurement 

The  power  deposited  at  a  specific  point  in  tissue  by  the  ultrasound  transducer  depends  on 
the  intensity  of  the  sound  field  at  that  point  and  the  absorption  coefficient  of  the  tissue  at  that 
point.  Clearly,  the  intensity  at  a  point  depends  on  the  initial  intensity  and  the  attenuation  of  the 
signal  in  the  tissue.  This  attenuation  is  due  to  scattering,  absorption,  and  geometric  attenuation. 
Previous  research  has  shown  that  the  attenuation  is  highly  variable  within  human  body  tissue,  with 
measurements  varying  significantly  for  different  tissue  types,  for  similar  tissue  on  different 
subjects,  and  even  for  the  same  tissue  on  the  same  patient  on  different  days.  [6-9]  Thus,  it  is 
important  to  measure  the  attenuation  of  the  ultrasonic  signal  in  vivo  on  the  patient  at  the 
treatment  site. 

Using  the  data  acquisition  system  shown  in  Fig.  2,  it  is  possible  to  estimate  the  attenuation 
of  the  ultrasound  at  different  points  in  the  patient  tissue  from  the  RF  signal  collected.  This  is 
accomplished  using  a  Log  Spectral  Difference  techmque.[7,8,10]  Details  of  this  technique  are 
discussed  in  the  second  year  progress  report.  It  is  important  to  note  that  Log  Spectral  Difference 
is  one  of  several  possible  techniques  and  provides  an  estimate  of  the  attenuation.  [5]  For  this 
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reason  it  is  important  to  test  the  accuracy  of  the  technique  with  other  measurements. 

Tests  were  conducted  where  the  acoustic  attenuation  of  a  series  of  samples  were 
measured  using  both  log  spectral  difference  and  the  Force  Balance  technique.  Using  the  Force 
Balance  technique  the  force  generated  when  an  acoustic  wave  impinges  on  an  absorbing  material 
is  measured  and  related  to  the  power  of  the  wave.  Using  this  technique  it  is  possible  to  make  a 
direct  measurement  of  the  attenuation  or  absorbed  power.  For  this  reason,  the  force  balance 
technique  is  used  as  a  standard  to  which  the  log  spectral  difference  measurements  are  compared. 
While  it  is  generally  considered  accurate  and  reliable,  the  force  balance  technique  is  not  suitable 
for  use  in  the  clinical  experiments  due  to  geometric  constraints. 

Using  the  techniques  described  above,  initial  tests  have  been  run  on  chicken  breast  to 
compare  the  results  of  the  two  techniques  discussed  above.  Using  the  Log  Spectral  Difference 
Technique,  the  normalized  attenuation  coefficient  for  the  chicken  was  measured  as  0.1608  Nepers 
x  cm 1  x MHz'1. [4]  Using  the  Force  Balance  system  measured  values  ranged  from  0. 108  to  0.138 
Nepers  x  cm1  x  MHz1.  These  values  are  of  the  same  magnitude  as  values  reported  in  the  literature 
(0.06  -  0. 13  Nepers x  cm'1  x MHz'1).^  1] 

Subsystem  2:  Absorbed  Power  Measurement 
No  work  accomplished  in  year  2. 

Subsystem  3:  Ultrasound  Propagation  Modeling 

During  the  second  year  two  competing  models  were  developed  and  compared.  The  first 
model  was  a  hybrid  numerical  model  using  integral  methods  to  solve  for  the  acoustic 
displacements  in  the  region  where  the  ultrasound  propagates  from  the  transducer  to  the  patient’s 
skin  and  using  a  finite  element  solution  of  the  wave  equation  for  propagation  in  the  body.  The 
model  predicts  the  displacements  and  stresses  propagating  through  the  tissue,  automatically 
accounting  for  reflections  and  transmission  at  all  interfaces.  This  model  was  expected  to  be  quite 
accurate.  A  second  model,  based  on  the  parabolic  equation  was  also  developed.  This  model 
accounts  for  forward  scattering  and  refraction,  but  ignores  reflections.  It  was  expected  to  be 
much  faster,  but  it  was  unclear  if  this  model  would  provide  sufficient  accuracy. 

Results:  Hybrid  Model 

The  Hybrid  Finite  Element  model  presented  two  difficulties: 

1 .  To  keep  the  computational  time  reasonable,  it  was  necessary  to  model  a  relatively 
small  volume.  It  was  not  possible  to  produce  boundary  conditions  at  the  edge  of  the 
volume  that  approximated  free  space  and  the  model  results  included  significant 
standing  waves  generated  at  the  model  boundaries.  Expanding  the  model  to  include 
the  entire  chest  cavity  would  solve  this  problem,  but  the  model  would  become 
unmanageably  large. 

2.  Initial  tests  were  performed  at  reduced  frequency  with  success.  When  frequencies 
corresponding  to  the  treatment  transducers  were  used,  the  finer  mesh  required  resulted 
in  extremely  high  computation  times. 
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Results:  Parabolic  Model 

The  parabolic  equation  method  is  an  efficient 
method  for  modeling  acoustic  wave  propagation  through 
low  contrast  acoustic  materials.  This  solution  technique  is 
somewhat  more  of  an  approximation  than  the  hybrid  model, 
but  because  of  the  geometry  used  in  the  hyperthermia 
experiment,  and  the  relatively  low  contrast  of  breast  tissue, 
it  is  possible  to  utilize  this  efficient  approximation.  The 
parabolic  model  was  found  to  be  computationally  8  time 
more  efficient  than  the  Green’s  function  model  and  64  times 
more  efficient  than  the  hybrid  model. 

Figures  3  and  4  demonstrate  the  kind  of  information 
available  from  the  models.  In  Figure  3  we  see  a  continuous 
plane  wave  propagating  from  the  left,  entering  the  body 
through  a  vertical  skin  layer.  The  wave  travels  through  an 
elliptical  fatty  region  and  into  the  spherical  tumor.  The  top 
figure  shows  the  wave  energy  density  that  tends  to  get 
weaker  as  the  wave  moves  to  the  right.  Note  that  the 
tumor  absorbs  most  of  the  wave  energy  resulting  in  a 
“dark”  shadow  behind  the  tumor.  The  bottom  figure  shows 
the  energy  absorbed.  Note  that  in  the  water  and  normal 
tissue  the  absorbed  energy  is  low;  even  when  the  wave 
energy  is  high.  There  is  an  increase  in  absorbed  energy  as 
the  wave  passes  into  the  tumor.  This  results  from  the  higher 
absorption  coefficient  in  the  tumor.  This  effect  is  shown 
more  clearly  in  Fig.  4  which  is  a  line  plot  of  the  absorbed 
energy  along  a  central  axis  in  Fig.  3. 

Model  Accuracy 

To  test  the  accuracy  of  the  parabolic  method,  results  were  compared  to  an  integral 
equation  solution  of  the  same  problem  that 
is  representative  of  the  best  models 
currently  in  use.  The  following  figures 
demonstrate  the  agreement  between  the 
integral  equation  method  and  the  parabolic 
method.  Figure  5  shows  the  geometry 
used  for  this  analysis.  Figure  6a  shows  a 
color  mapped  representation  of  a  slice  of 
the  field.  Figure  6b  is  an  axial  plot  of  the 
same  data.  Figure  7a  and  7b  show  the 
results  using  a  speed  of  sound  in  the  sphere 
significantly  less  than  and  significantly 
greater  than  that  of  water.  Given  the 
differences  in  the  two  models,  one  would 
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Figure  3:  The  upper  image  is  the  incident 
energy  density.  The  lower  image 
shows  energy  absorbed  by  the  tissue. 
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expect  error  to  increase  as  the  mismatch  in  sound 
velocity  between  the  tumor  and  water  increases. 

Thus,  these  two  conditions  should  represent  worst 
case  errors. 

The  excellent  match  between  the  two  models 
in  the  forward  scattered  diffraction  patterns  shown 
in  Fig.  7  and  the  poor  match  in  the  reflected 
diffraction  pattern  shown  in  Fig.  6b  highlight  the 
differences  between  the  two  models.  The  excellent 
overall  match  between  the  two  models  in  these 
worst  case  conditions,  coupled  with  the  speed 
provided  by  the  parabolic  model  make  it  an  excellent  candidate  for  our  application. 


Comparison  of  axial  field 


Figure  6:  a)  Color  mapped  image  of  the  power  deposition  in  the  model  from  Fig.  5. 

b)  The  integral  equation  is  shown  in  solid,  and  shows  the  diffraction  pattern  due  to  reflections  from  the 
hard  sphere.  This  plot  shows  quantitatively  that  the  parabolic  method  does  not  take  into  account  the 
reflections  of  sound  waves.  Conditions  for  this  graph  correspond  to  those  in  Fig.  7b. 


Figure  7:  a)  The  speed  of  sound  internal  to  the  sphere  was  1267.7  m/sec,  substantially  less  than  that  which  would 
be  encountered  in  the  breast,  thereby  making  the  impedance  contrast  quite  a  bit  larger  than  is  likely  to  be 
encountered  in  practice,  b)  Here,  the  speed  of  sound  internal  to  the  sphere  was  1677  m/sec,  which  is 
substantially  greater  than  any  tumor  that  would  be  found  in  the  breast.  In  each  case  the  solid  line  represents 
the  integral  model  and  the  dashed  line  represents  the  parabolic  model. 
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Subsystem  4:  Inverse  Optimization  Technique 
No  work  accomplished  during  the  2nd  year. 

Third  Year  Accomplishments 

Subsystem  1:  Ultrasound  Based  Anatomy  and  Attenuation  Measurement 

•  Laboratory  hardware  was  developed  to  provide  3D  scanning  capability  for  the  medical 
B-scan  imager.  Clinical  hardware  is  now  essentially  complete. 

•  Software  to  produce  a  3D,  tissue  attenuation  map  from  the  A-mode  signals  was 
completed  and  is  undergoing  verification  testing. 

•  Software  and  hardware  for  clinical  patient  geometry  acquisition  is  approximately  50% 
complete. 

Patient  Geometry  Acquisition  and  Registration 

The  goal  here  is  to  be  able  to  place  a  patient  in  a  known,  repeatable  location  on  the 
treatment  station  and,  using  the  B-scan  imager,  identify  the  geometric  structure  of  the  tissue  and 
tumor  in  the  region  under  treatment.  The  proposed  solution  to  this  problem  has  evolved 
significantly  since  the  initial  proposal.  In  the  initial  proposal  CT  data  and  printed  B-scan  images 
were  to  be  converted  into  a  geometric  description  of  the  treatment  region  suitable  for  use  by  the 
math  model.  This  has  evolved  into  specific  hardware  and  software  that  make  up  a  clinical  system 
that  is  described  below.  The  hardware  for  this  system  was  not  included  in  the  original  proposal 
and  we  are  currently  using  borrowed  equipment  and  awaiting  approval  to  purchase  equipment 
that  will  be  dedicated  to 
this  work.  During  this 
reporting  period,  the 
software  for  the  system 
described  below  was 
developed  and  is 
currently  approximately 
50%  finished. 


Clinical  System 
Description 

The  patient 
positioning  system 
consists  of  several 
integrated  steps.  The  first 
step  is  to  reproduce  the 
laser  positioning  system 
used  for  radiation 
therapy.  Use  if  this 
positioning  system  will 
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Figure  8:  The  software  under  development  acquires  a  series  of  B-scan  images, 
tiles  them  together  creating  a  larger  2D  image,  and  coordinates  the 
display  of  the  2D  images  to  represent  3D  geometry. 
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allow  CT  scans  taken  by  the  radiologist  to  be  used  for  added  reference.  The  system  requires  three 
collimated,  wall-mounted  lasers  that  are  aligned  with  points  marked  by  freckle-sized  tattoos  on 
the  patient’s  skin.  After  the  patient  is  initially  positioned,  the  laser  spots  on  the  patient’s  skin  are 
marked  with  tattoos.  The  patient  position  is  checked  in  later  treatments  by  observing  the  position 
of  the  laser  spots  with  respect  to  the  tattoos.  The  hyperthermia  treatment  apparatus  at  the 
University  of  Utah  currently  is  equipped  with  two  wall-mounted  lasers  and  a  third  is  planned. 

In  the  second  phase  of  positioning  the  patient,  a  series  of  B-scan  ultrasound  images  of  the 
treatment  area  are  collected.  These  images  are  used  to  identify  the  locations  of  thermocouples, 
interfaces  between  different  tissue  types,  and  other  anatomical  features  within  the  region  that  will 
be  affected  by  the  hyperthermia  treatment.  The  B-scan  images  are  digitized  using  a  video  capture 
card.  A  custom  software,  developed  as  part  of  this  research,  coordinates  multiple  B-scan  images 
into  a  series  of  slices  representing  the  3D  patient  geometry.  A  screen  capture  from  this  system  is 
shown  in  Figure  8. 

The  doctor  or  technician  traces  features  of  the  B-scan  images  that  are  of  interest.  Using 
an  attenuation  map  developed  using  a  process  described  below  as  a  reference,  the  different 
regions  identified  by  the  tracings  are  assigned  attenuation  values.  The  software  processes  this 
information  and  produces  a  3D  grid  of  nodes  where  each  node  is  assigned  appropriate  values  to 
represent  the  patient  geometry  for  processing  by  the  math  model.  Figure  9  demonstrates  this 
process. 

The  left  most  image  contains  3  B-scans  that  have  been  registered  and  combined, 
producing  a  single  image.  In  the  center  image,  the  doctor  or  technician  has  mouse  sketched 
different  regions.  In  the  right  most  image,  these  regions  have  been  assigned  different  acoustic 
properties  based  on  apriori  knowledge  and  information  from  the  attenuation  mapping  system 
discussed  below. 

The  following  outline  describes  the  anticipated  clinical  procedure  using  the  developed 
system.  Depending  on  the  time  required  to  complete  the  inverse  optimization  procedure,  the 
patient  may  either  lay  in  position  and  wait,  or  be  released  and  return  for  treatment  hours  or  days 
later. 


Figure  9:  Each  image  is  processed  by  a  doctor  or  technician.  The  process  involves  mouse  sketching  the  borders  of 
identifiable  regions  and  assigning  each  region  uniform  acoustic  values  estimated  from  an  attenuation  map. 
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Patient  geometry  acquisition  and  initial  registration: 

1 .  A  coordinate  system  with  a  zero  point  near  the  tumor  is  usually  defined  by  the  radiologist. 
This  zero  point  can  be  located  by  aligning  wall-mounted  lasers  with  three  small  tattooed 
marks  on  the  patients  skin  or  with  marks  on  a  mold  used  to  hold  the  patient  in  position  during 
the  treatment.  Align  the  patient  at  the  zero  point. 

2.  Center  the  B-scan  around  this  zero  point  and  scan  the  transducer  in  order  to  capture  a  “grid” 
of  two-dimensional  images  that  have  known  locations  with  respect  to  the  zero  point. 

3 .  Use  a  frame  grabber  to  store  the  images.  The  location  of  the  scan  with  respect  to  the  zero 
point  must  be  recorded  with  each  image. 

4.  Trace  interfaces  and  loops  that  represent  different  tissue  regions.  Mark  thermocouple 
positions. 

5.  Store  these  traces  in  an  image  that  has  the  same  scale  and  position  as  the  original  image. 

6.  Use  attenuation  information  to  define  the  tissue  properties  within  the  traced  regions. 

7.  Format  data  and  output  to  optimization  software. 

Repositioning  patient  for  treatments: 

1 .  Align  wall-mounted  lasers  with  marks  on  the  patient’s  skin. 

2.  Scan  the  anatomical  region  to  be  treated  with  an  ultrasound  B-scan  transducer  using  the  same 
grid  pattern  used  in  the  initial  geometry  acquisition  scan. 

3.  Use  a  frame  grabber  to  store  the  images  and  their  location  with  respect  to  the  zero  point. 

4.  Compare  the  captured  images  to  those  from  previous  scans  to  verify  that 
thermocouples  and  markers  on  the  skin  are  in  the  same  locations  in  the  corresponding  images. 
Adjust  patient  position  until  the  images  show  that  the  thermocouples  and  markers  are  in  the 
correct  locations. 
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Three-dimensional  Ultrasonic  Attenuation  Field  Measurement 


As  mentioned  above,  the  goal  of  this  research  is  to  produce  a  model  that  will  accurately 
predict  the  power  deposited  at  all  points  in  a  region  of  tissue,  during  a  hyperthermia  treatment. 
This  will  allow  the  physician  to  develop  a  better  treatment  plan  that  will  maximize  the  positive 
effects  and  minimize  pain  and  other  adverse  effects.  In  order  for  the  inverse  model  to  accomplish 
this  efficiently,  the  model  must  begin  with  model  parameters  that  are  “in  the  ball  park”  of  the  final 
optimal  values.  The  least  well  known  of  these  parameters  is  the  acoustic  attenuation  of  the  tissue. 
Thus,  a  technique  is  developed  here  that  will  predict  the  acoustic  attenuation  field  within  the 
tissue  based  on  data  from  a  B-scan  imager.  The  third  year  efforts  have  extended  the  previous 
accomplishments  to  include  the  ability  to  produce  a  3D  attenuation  map  from  a  series  of  B-scan 
images. 

Background 

Measurement  of  ultrasonic  attenuation  from  reflected  signals  has  shown  promise  of  being 
a  valuable  diagnostic  tool[12-20].  It  allows  physicians  to  diagnose  illnesses  that  cannot  be 
diagnosed  on  a  B-mode  ultrasound  image.  Because  of  this,  several  methods  of  calculating  the 
ultrasonic  attenuation  from  reflected  signals  have  been  developed  [14,17,19-22],  These  methods 
were  then  used  to  calculate  the  global  attenuation  of  a  structure.  The  next  step  in  attenuation 
calculation  was  to  produce  an  attenuation  image  that  could  be  used  in  conjunction  with  the  B- 
mode  image.  Walach  and  others  have  developed  techniques  that  relate  a  two-dimensional  B- 
mode  image  to  a  two-dimensional  attenuation  image  or  map  [12,13,23],  Seeing  both  images  side- 
by-side,  the  physician  has  an  increased  diagnostic  capability. 

However,  the  usefulness  of  an  attenuation  map  is  not  limited  to  diagnostics.  A  3D 
attenuation  mapping  technique  has  been  developed  here  that  not  only  provides  the  initial  estimates 
for  the  power  deposition  model.  But  has  potential  to  produce  enhanced  resolution  2D  attenuation 
images. 

Method 

In  order  to  build  the  three-dimensional  attenuation  map,  a  B-mode  imaging  transducer  was 
placed  on  a  two-dimensional  linear  motion  machine.  The  transducer  is  part  of  an  ultrasonic 
imaging  system  (PIE  Scanner  250)  which  provides  access  to  the  reflected  ultrasound  signal  (A- 
line),  before  processing,  as  well  as  synchronization  signals  to  give  the  location  of  the  waveform 
along  the  width  of  the  transducer.  With  this  in  place,  software  was  developed  that  scanned  the 
transducer  to  cover  the  entire  geometry  of  the  structure  being  modeled,  while  collecting  and 
organizing  the  data.  This  is  the  same  data  shown  in  Figures  8  and  9  collected  by  the  patient 
registration  system,  except  that  the  raw  A-lines  are  recorded. 

The  data  was  collected  with  an  eight-bit  analog  to  digital  converter,  and  was  then  written 
to  disk.  A  technique  has  been  developed  that  uses  several  different  gain  settings  for  each 
collection  to  mimic  a  higher-precision  conversion.  Using  this  technique,  each  A-line  was  analyzed 
at  the  highest  possible  gain  setting  where  saturation  did  not  occur.  The  focus  of  the  transducer 
was  set  so  that  the  structure  being  analyzed  was  in  the  far  field  of  the  transducer.  This  insures 
that  diffraction  effects  in  the  near  field  do  not  enter  into  the  calculations.  Once  the  data  was 
collected  it  was  then  analyzed  to  produce  the  attenuation  map. 
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The  log  spectral  difference 
method  was  used  to  estimate  the  slope 
of  the  attenuation  to  frequency 
relationship.  (Lyons  and  Parker  found 
that  for  most  soft  tissues,  the 
attenuation  was  related  to  frequency 
to  the  power  of  n,  where  n  is  equal  to 
1 .0  - 1 .3  [25],  The  value  for  n  was 
set  to  1  in  this  research.)  This  method 
utilizes  the  ratio  of  the  power 
spectrum  from  a  shallow  and  deep 
segment  of  the  region  of  interest  (see 
Figure  10)  to  compute  the  attenuation 
slope  for  that  region.  The  specific 
steps  of  this  technique  have  been 
presented  in  the  2nd  year  progress 
report  as  well  as  previously  in  the 
literature  and  will  not  be  discussed 
here  [14,15,22,26],  Our  application  of  this  method  for  the  production  of  a  3D,  attenuation  map 
will  be  discussed  next. 

As  mentioned  previously,  the  log  spectral  difference  technique  is  a  way  of  estimating  the 
slope  of  the  attenuation  for  a  specified  Region  Of  Interest  (ROI).  To  produce  the  attenuation 
map,  a  ROI  of  a  determined  size  was  moved  throughout  the  three-dimensional  data  space  until  an 
attenuation  map  was  produced.  The  region  of  interest  selected  was  1 .0  cm  wide  and  1 .7  cm  in 
depth  (see  Figure  10).  The  width  was  chosen  because  it  allowed  the  computation  to  be  averaged 
over  ten  reflected  signals.  The  depth  for  the  region  of  interest  came  from  a  formula  presented  by 
Kuc  [26]  and  the  specifications 
of  the  transducer.  Once  the  ROI 
size  was  selected,  it  was  moved 
throughout  the  data  space 
incrementally  in  1.0  mm  steps  in 
all  three  directions.  The 
attenuation-slope  value  for  each 
calculation  was  assigned  to  a 
point  in  the  center  of  the  ROI, 
thus  producing  an  attenuation 
map.  This  2D  ROI  is  similar  to 
that  used  by  previous 
researchers.  The  additional 
information  available  from  the 
3D  data  collection  used  here 
allows  a  3D  ROI  to  be 
investigated,  where  data  from  a 
3D  region  are  averaged  together 
to  improve  the  accuracy  and/or 
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resolution  of  the  measurement.  Software  has  been  developed  to  completely  automate  the  task  of 
data  acquisition  and  analysis,  and  calibration  tests  are  underway. 


Results 

Various  tests  have  been  performed  to  determine  the  settings  that  produced  the  most 
precise  and  least  noisy  attenuation  map.  Initial  tests  varied  the  ultrasound  power,  gain  setting, 
and  whether  or  not  the  region  of  interest  should  be  two-dimensional  (space  averaged)  or  three- 
dimensional  (volume  averaged).  For  these  preliminary  experiments,  visual  comparison  of  the 
attenuation  images  versus  the  anticipated  attenuation  image  (obtained  from  examining  the  B-mode 
image)  was  used  to  determine  the  optimal  settings. 

Preliminary  results  showed  no  benefit  from  the  3D  volume  averaging.  We  anticipate  that 
this  is  due  to  registration  problems  in  the  data  collection  and  we  are  continuing  to  investigate. 
During  the  analysis  of  the  attenuation  images,  it  was  noted  that  data  was  lost  in  areas  where  the 
ultrasound  back  scatter  was  too  low,  or  too  high.  For  any  single  gain  setting,  the  data  loss  was 
significant.  However,  based  on  a  visual  comparison,  other  gain  settings  appeared  to  produce 
attenuation  images  of  the  almost  the  same  quality,  but  with  data  losses  in  different  regions.  These 
observations  resulted  in  the  multiple  gain  setting  technique  discussed  above.  This  technique  has 
resulted  in  images  with  significantly  less  data  loss. 

Figure  1 1  presents  a  B-mode  image  of  a  breast  biopsy  phantom.  This  image  was  taken  at 
the  same  time  that  attenuation  data  was  being  acquired  for  the  same  image  space.  Because  the 
attenuation  slope  value  is  assigned  to  a  point  in  the  center  of  the  ROI,  and  the  because  of  the  size 
of  the  ROI,  0.5  cm  of  attenuation  data  was  lost  on  each  side  of  the  B-mode  image.  For  the  same 
reason,  0.85  cm  was  lost  from  the  top  and  bottom  of  the 
image.  Figure  12  shows  the  corresponding  attenuation  map 
for  the  region  indicated  in  Figure  1 1 . 

Discussion 

Attenuation  maps  are  not  a  new  application. 

Walach  and  others  have  produced  attenuation  maps  in  the 
past  [12,13,23],  However,  this  application  is  the  first 
three-dimensional  mapping  system  we  have  been  able  to 
find.  This  application  has  limitations  similar  to  the  other 
maps  that  have  been  produced.  One  such  limitation  is  the 
resolution  of  the  attenuation  map. 

Most  attenuation-calculation  techniques  require  a 
large  ROI  in  the  depth  direction.  The  literature  has 
recommended  that  the  log  spectral  difference  method  be 
used  for  global  attenuation  estimation  (greater  than  5.0  cm) 

[14,22],  However,  Kuc  has  theorized  that  this  method  can 
be  used  for  smaller  regions,  but  that  the  size  of  the  ROI  has 
a  lower  limit.  As  the  length  of  the  ROI  decreases  below  this 
limit,  the  noise  associated  with  the  application  increases 
dramatically  [26],  As  can  be  seen  from  Figure  12,  the 
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Figure  12:  Attenuation  image.  This 

image  corresponds  to  the  highlighted 
region  in  Figure  2.  Values  of 
attenuation  are  in  l/(cm  dB  MHz). 
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resolution  is  poor  compared  to  the  B-scan  image.  In  addition,  the  attenuation  image  only  partially 
matches  the  expected  image  (from  examining  the  B-mode  image  in  Figure  11).  Experiments  using 
the  radiation  force  balance  apparatus  developed  in  the  first  and  second  years  of  the  project,  are 
currently  under  way  to  quantify  the  error  between  the  attenuation  predicted  by  this  technique  and 
the  true  attenuation  of  a  sample. 

Inhomogeneous  regions  are  another  problem  with  this  technique.  Several  interfaces  and 
varying  structures  were  present  in  the  phantom  that  was  used  for  calculation.  The  development  of 
this  and  other  estimation  techniques  were  based  on  the  tissue  being  homogeneous  within  the  ROI 
[14,17,18,21,22],  This  problem  was  initially  ignored  in  this  study  in  order  to  discover  how  the 
method  would  handle  interfaces  and  different  tissues  within  the  ROI..  The  bright  circle  in  Fig.  1 1 
is  a  water  filled  sphere,  1 .2  cm  diameter.  Note  that  while  this  sphere  is  smaller  than  the  ROI  it 
does  show  up.  This  is  taken  as  a  sign  that  this  method  may  be  used  to  calculate  the  attenuation  in 
small  tissue  structures. 

Future  Efforts 

Clearly  interfaces  present  the  opportunity  for  a  variation  of  attenuation  in  the  two  regions 
separated  by  the  interface.  When  the  ROI  overlaps  the  interface,  the  attenuation  measurement 
can  be  expected  to  be  inaccurate.  (While  it  may  accurately  represent  the  average  of  the 
attenuation  within  the  ROI,  we  anticipate  using  the  attenuation  map  to  allow  a  technician  to  select 
a  representative  value  for  an  identified  region.  In  this  case  we  prefer  attenuation  values  that  are 
representative  of  a  single  region.)  We  are  currently  developing  an  algorithm  to  separate  “good” 
data  from  noisy  measurements  and  anticipate  that  this  will  produce  larger  regions  of  “lost”  data 
with  smaller  regions,  centered  within  homogeneous  regions,  where  the  results  are  relatively 
uniform.  This  technique  is  based  on  an  examination  of  a  histogram  of  attenuation  values  within  a 
defined  neighborhood  of  each  other.  Using  this  technique,  attenuation  values  are  only  considered 
valid  if  the  histogram  has  a  single  mode,  with  a  standard  deviation  less  than  a  threshold  value.  At 
interface  regions  bi-modal  histograms  or  large  standard  deviations  are  expected. 

Another  source  of  error  in  the  three-dimensional  map  is  the  location  error  of  the 
transducer.  In  tests,  the  largest  error  in  positioning  the  transducer  found  was  0.04  mm.  The 
beam  width  of  this  transducer  at  the  focal  point  is  1 .4  mm.  Therefore,  the  positioning  error  is  not 
considered  a  significant  source  of  error. 

Regardless  of  the  errors  involved,  the  attenuation  estimation  method  should  be  successful 
if  it  can  appropriately  predict  the  relationship  between  the  attenuation  and  the  various  structures 
within  the  phantom  or  anatomy,  (i.e.  which  regions  have  “higher”  attenuation  and  which  regions 
have  lower  attenuation.)  This  kind  of  an  approach  can  be  taken  because  of  the  iterative  nature  of 
the  inverse  model  that  will  use  these  values. 


Subsystem  2:  Experimental  SAR  Measurement 

•  Thermocouple  based  beam  plots  have  been  performed  to  calibrate  the  ultrasound 
power  field  in  water  for  specific  transducers. 

•  Experiments  have  been  conducted  in  water  and  will  begin  shortly  using  cast  agar 
phantoms. 

•  Preliminary  in-vivo  experiments  were  conducted  on  a  dog. 
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•  A  finite  difference  model  has  been  constructed  to  help  interpret  the  experimental 
measurements 

To  accomplish  the  inverse  optimization  of  the  theoretical  model,  the  results  of  the  model 
must  be  compared  against  actual  values  that  are  determined  by  experimental  means.  The 
ultrasound  propagation  model  predicts  the  incident  power  at  all  points  in  the  modeled  region,  as 
well  as  the  absorbed  power.  At  a  given  point,  these  two  quantities  are  related  through  the  local 
attenuation.  Thus,  the  model  can  be  optimized  by  comparing  the  propagation  model  results  to 
either  the  incident  power  or  the  absorbed  power  at  several  points  in  the  model.  The  experimental 
values  needed  for  optimization  can  be  determined  using  the  same  thermocouples  which  are 
inserted  for  observation  of  steady  state  temperature  during  clinical  hyperthermia  treatment. 

Experimental  measurement  of  the  power  deposited  in  tissue  at  a  specific  point  is  not 
entirely  straight  forward.  Researchers  note  that  the  interaction  of  the  thermocouple  probe  with 
the  ultrasound  results  in  temperature  artifact.  [24]  This  artifact  results  from  two  sources:  viscous 
heating  at  the  interface  between  the  probe  surface  and  the  tissue,  and  absorption  of  acoustic 
energy  by  the  probe  its  self.  The  viscous  heating  is  a  surface  phenomena  and  results  in  a  small 
artifact  with  a  short  time  constant  causing  a  rapid  initial  temperature  rise  when  acoustic  power  is 
initiated.  This  artifact  is  usually  only  noticed  when  conducting  tests  using  bare  thermocouples. 
Probe  sheaths  of  steel,  fused  silica,  and  various  plastics  have  been  shown  to  absorb  ultrasound 
power  at  a  higher  rate  than  tissue.  Plastics,  the  most  common  sheath  materials  used  in 
hyperthermia  treatment,  absorb  ultrasound  energy  at  a  significantly  higher  rate  than  tissue 
resulting  in  a  significant  artifact  that  is  exhibited  over  a  longer  time  scale  than  the  viscous  artifact. 
The  temperature  rise  in  the  tissue,  a  much  larger  volume  than  the  sheath,  occurs  over  an  even 
longer  time  scale. 

The  combined  result  of 
these  effects  is  shown  in  Figure 
13.  In  Figure  1 3  (a),  a  bare 
thermocouple  was  submersed  in 
a  water  bath  and  insonified  for 
0.6  seconds.  Note  the  rapid 
initial  jump  followed  by  a 
relatively  linear  region.  The 
initial  jump  is  caused  by  viscous 
heating  at  the 

thermocouple/water  interface 
and  the  linear  region  is 
established  when  local 
conduction  between  the 
thermocouple  and  surrounding 
water  becomes  significant.  At 
longer  times  the  temperature 
rise  will  become  non-linear  as 
larger  scale  conduction  and 
convection  become  significant. 

In  Figure  13(b)  the  experiment 


BARE  THERMOCOUPLE  SHEATHED  THERMOCOUPLE 
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Figure  13 :  a)  Results  from  a  bare  thermocouple  in  water  insonified  for  0.6 
seconds,  b)  Results  from  a  plastic  sheathed  thermocouple  in 
water  insonified  for  0.6  seconds. 
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was  repeated  with  a 
sheathed  thermocouple. 

Note  that  the  temperature 
rise  is  much  more 
significant,  rising  roughly  20 
degrees  in  0.6  seconds 
compared  to  less  than  one 
degree  for  the  bare 
thermocouple.  This  large 
jump  is  attributed  to  the 
power  absorption  of  the 
sheath.  Any  viscous  heating 
that  occurs  is  clearly  lost  in 
this  rapid  temperature  rise. 

The  slope  of  this  rise  is 
indicative  of  the  power 
absorption  in  the  sheath. 

If  the  power  is 
applied  for  a  longer  time, 
the  rapid  initial  rise  of  the 
sheathed  thermocouple  is 
balanced  by  conduction  between  the  sheath  and  the  surrounding  water.  This  effect  is  shown  in 
Figure  14.  The  second  linear  regime  achieved  represents  the  temperature  rise  of  the  media 
surrounding  the  thermocouple. 

In  these  experiments  the  thermocouples  were  immersed  in  water,  a  very  poor  absorber  of 
ultrasound  energy.  Thus,  essentially  all  of  the  temperature  rise  of  the  water  was  due  to 
conduction  from  the  thermocouple  sheath.  In  tissue  we  expect  larger  absorption  directly  by  the 
tissue.  An  actual  experiment  conducted  in  tissue  will  produce  results  similar  to  those  shown  in 
Figure  15.  Here  a  sheathed  thermocouple  was  inserted  in  the  muscle  tissue  of  the  thigh  of  a  dog. 
The  focus  of  the  ultrasound  transducer  was  centered  on  the  thermocouple  and  the  area  was 
insonified  for  0.6  seconds  at  several  different  power  levels.  Note  that  in  each  case  the  power  was 
much  lower  than  that  used  in  the  experiment  from  Figure  14,  evidenced  by  the  relatively  low 
temperature  rise.  In  future  experiments  we  intend  to  calibrate  the  plastic  sheath  and  use  it  as  an 
incident  power  sensor,  allowing  the  ultrasound  power  incident  at  the  thermocouple  to  be 
measured. 


Modeling  the  Thermocouple  Behavior 


In  the  absence  of  conduction,  power  is  related  to  heat  generation  by  the  following 
equation: 


dT 

^T,=cd 


Equation  1 


where  p  is  the  density,  Cp  is  the  specific  heat,  a  is  the  absorption  coefficient,  and  I  is  the  local 
intensity  of  the  ultrasound  field. 
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rate  of  heating  is  0.68  C/s  . 


rate  of  heating  is  0,77  C/s 
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Impulse  test  of  dog  at  "35%"  power 


Impulse  test  of  dog  at  "40%"  power 


Impulse  test  of  dog  at  "45%"  power 
rate  of  heating  is  1 .03  C/s 

_ _ 5 

Impulse  test  of  dog  at  "50%"  power 
rate  of  heating  is  1.39  C/s 
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Figure  15:  Four  sequential  tests  of  a  plastic  sheathed  thermocouple  inserted  in  the  thigh  of  a 
dog.  _ 
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If  a  plastic  thermocouple  probe  and  surrounding  tissue  at  uniform  temperature  is  suddenly 
insonated  by  ultrasound  of  constant  intensity,  the  probe  and  tissue  will  begin  to  rise  in 
temperature.  The  plastic  probe  will  heat  more  rapidly  than  the  tissue  because  it  has  a  larger 
absorption  coefficient  than  the  tissue.  For  a  short  period  of  time  immediately  following 
insonation,  conduction  of  heat  will  be  small  and  the  rate  of  temperature  rise  is  dependent  solely  on 
the  absorption  coefficient  of  the  medium  being  insonated.  If  the  time  constant  of  the 
thermocouple  is  sufficiently  small  the  thermocouple,  being  surrounded  by  plastic,  will  reflect  the 
temperature  of  the  plastic.  If  the  absorption  characteristics  of  the  thermocouple  probe  can  be 
sufficiently  characterized,  the  intensity  of  incident  ultrasound  can  be  determined  by  the 
temperature  vs.  time  response  of  the  thermocouple  to  a  suddenly  applied  ultrasound  field. 

The  radial  distribution  of  a  focussed  ultrasound  beam  can  be  closely  estimated  as  a 
Gaussian  function.  The  following  equation  describes  the  radial  intensity  distribution: 

I(r)  =  Im y'"  Equation  2 

where  Imax  is  the  intensity  at  r  =  0  and  P  is  a  measure  of  the  beam  width  as  determined  by  a  least 
squares  error  curve  fit.  [24] 

If  the  ultrasound  is  applied  for  a  short  duration  to  prevent  significant  conduction  effects 
then  the  temperature  distribution  produced  by  the  absorption  of  the  ultrasound  will  also  be 
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Gaussian.  Using  the  theory  put  forth  by  Parker,  the  time  dependent  temperature  at  the  center  of 
the  focal  region  can  be  described  by  the  following  equation: 

T(t)  =  Tmax  l(4kt  1/3  +  1)  Equation  3 

where  k  is  the  heat  conduction  coefficient  of  the  tissue  and  7max  is  the  maximum  temperature 
attained  at  the  center  of  the  focal  region  at  the  end  of  ultrasound  application. 

At  the  termination  of  the  momentary  application  of  ultrasound  the  plastic  sheathing  of  the 
thermocouple  probe  will  be  at  a  higher  temperature  than  the  surrounding  tissue.  However, 
because  the  heat  capacity  of  the  thermocouple  is  small  compared  to  the  heat  capacity  of  the  tissue, 
excess  heat  stored  in  the  probe  will  diffuse  into  the  tissue  and  the  probe  temperature  will  reach 
equilibrium  with  the  surrounding  tissue  and  follow  the  temperature  decay  of  the  surrounding 
tissue.  Curve  fitting  the  decay  of  the  tissue  temperature  by  the  least  squares  error  method  the 
decay  can  be  backtracked  and  the  value  of  Tmax  determined.  Tmax  and  the  duration  of  insonation 
are  used  with  equation  (1)  and  the  intensity  determined  by  the  rate  of  heating  to  determine  the 
absorption  coefficient  of  the  tissue.  [24] 

It  has  been  determined  by  previous  investigators  that  absorption  accounts  for  90-100%  of 
attenuation  in  soft  tissue.  [25]  Therefore  it  is  not  completely  accurate  to  compare  attenuation 
determined  by  B-mode  imager  and  absorption  determined  by  thermocouple  probe  but  an  error  less 
than  10%  is  not  unreasonable  when  compared  to  other  errors  likely  to  exist. 

A  program  was  created  on  the  hospital  treatment  system  that  positions  the  transducer, 
delivers  power  for  a  period  of  time,  and  records  the  temporal  temperature  data  from  the 
thermocouple  probes.  This  program  is  able  to  attain  temperature  data  from  several  seven  junction 
probes  every  3/100  sec.  To  test  the  program  a  10  cm,  2MHz  transducer  was  mounted  to  the 
positioning  gantry  and  an  array  of  thermocouples  was  placed  in  a  water  bath  above  the  treatment 
system.  The  thermocouples  were  located  with  the  existing  thermocouple  locate  routine.  It  was 
found  that  the  positioning  system  was  only  able  to  locate  the  thermocouples  with  ±  2-4cm 
accuracy.  A  radial  beam  intensity  plot  was  attempted  beginning  1  cm  away  from  the  center  of 
focus  and  traversing  through  the  focus  at  2mm  increments.  Because  of  the  tight  beam  focus  and 
lack  of  positioning  accuracy  it  was  found  that  there  was  very  little  correlation  between  expected 
position  and  beam  intensity. 

To  characterize  the  beam  shapes  of  the  available  transducers  an  array  of  epoxy  covered 
thermocouples  was  created  and  oriented  perpendicular  to  the  beam  axis.  The  array  was  moved 
through  the  beam  transverse  to  the  beam  axis  in  1mm  increments  and  the  steady  state  temperature 
was  recorded  for  the  thermocouples  at  several  different  distances  from  the  face  of  the  transducer. 
This  procedure  was  followed  for  the  focussed  10cm  diameter  2MHz,  focussed  7cm  diameter 
500kHz,  unfocussed  square  4cm  1MHz,  and  unfocused  3  cm  diameter  4MHz  transducers.  The 
unfocussed  transducers  were  not  characterizable  because  of  the  irregular  intensity  consistent  with 
the  long  near  field  of  the  beams.  The  10cm  and  7cm  diameter  transducers  had  a  Gaussian  profile 
in  the  radial  direction  at  the  focus.  Least  squares  error  curve  fit  produced  a  Gaussian  parameter 
of  0.029  cm'2  for  the  10  cm  diameter  transducer  and  a  Gaussian  parameter  of  .42  cm'2  for  the  7cm 
diameter  transducer. 
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The  0.029  cm'2  Gaussian  parameter  at  the  focus  of  the  10cm  diameter  transducer  creates  a 
beam  width  of  around  6mm.  which  is  too  narrow  for  the  accuracy  of  the  hospital  positioning 
system  and  the  steep  temperature  gradients  created  by  the  narrow  focus  decrease  the  temporal 
window  during  which  conduction  is  insignificant.  The  0.42  cm'2  Gaussian  parameter  at  the  focus 
of  the  7  cm.  diameter  transducer  creates  a  beam  width  larger  than  2  cm.  This  larger  beam  width 
will  create  a  less  steep  temperature  gradient  and  lengthen  the  temporal  window  for  the  rate  of 
heating  measurement.  The  longer  length  of  focus  of  the  7  cm.  transducer  is  also  beneficial  in 
avoiding  error  in  the  Z-direction. 

Finite  difference  code  was  written  to  model  a  thermocouple  probe  embedded  in  tissue. 
The  finite  difference  model  will  be  used  to  verify  the  assumptions  and  clarify  the  unknowns  of  the 
experimental  technique.  Verification  of  the  1-D  model  is  in  process  after  which  a  2-D  and/or  3-D 
model  will  be  created. 

A  considerable  amount  of  phantom  testing  with  the  7cm.  diameter  transducer  will  begin 
shortly.  The  phantom  tests  will  be  designed  to  prove  the  accuracy  of  the  rate  of  heating  and 
thermal  pulse  decay  technique  in  combination  with  the  probe  geometry  used  in  clinical 
hyperthermia  treatment.  After  verification  the  procedure  will  be  integrated  with  the  other 
subsystems  of  the  inverse  estimation  technique  to  provide  geometric  data  for  the  thermocouple 
probes  and  local  measured  power  and  absorption  values. 

Subsystem  3:  Model  Development 

•  The  parabolic  model  has  been  selected  for  use  in  regions  of  low  acoustic  contrast.  In 
these  regions  the  model  is  quite  fast  and  loses  little  accuracy  when  compared  to  the 
Green’s  Function  Model.  The  Green’s  Function  Model  can  be  used  in  high  contrast 
areas  (close  to  ribs,  etc.) 

•  An  interface  has  been  built  for  the  model  which  allows  the  user  to  specify  the  volume 
modeled,  the  acoustic  power  input  as  a  2D  complex  pressure  function,  and  the  model 
parameters;  density,  sound  velocity,  and  absorption  at  each  node  within  the  volume. 

Based  on  the  results  of  the  2nd  year,  the  parabolic  model  has  been  selected  for  use  in 
regions  of  low  acoustic  contrast.  In  these  regions  the  model  is  quite  fast  and  loses  little  accuracy 
when  compared  to  the  Green’s  Function  Model.  If  necessary,  the  Green’s  Function  Model  can 
be  used  in  high  contrast  areas  (close  to  ribs,  etc.)  The  Green’s  Function  Model  could  be  used 
everywhere  resulting  in  a  significantly  slower,  slightly  more  accurate  result. 

During  the  third,  year  work  on  the  model  has  focused  on  the  development  of  a  practical 
and  convenient  user  interface  that  will  allow  the  model  to  be 
adapted  to  different  transducers  and  tissue  geometries.  Toward 
this  purpose,  an  interface  has  been  built  for  the  model  which 
allows  the  user  to  specify  the  volume  modeled,  the  acoustic 
power  input  as  a  2D  complex  pressure  function,  and  the  model 
parameters,  density,  sound  velocity,  and  absorption  at  each 
node  within  the  volume.  The  model  interface,  as  described 
below  is  complete. 
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The  interface  allows  the  user  to  specify  the  dimensions,  X,  Y,  Z  of  a  rectangular  volume 
that  will  be  modeled.  The  size  of  the  volume  to  be  modeled  is  limited  by  the  requirement  for  2 
nodes  per  wavelength  at  the  transducer  frequency  under  investigation  and  by  the  maximum 
number  of  nodes  that  the  model  can  handle.  At  the  current  time,  the  model  can  handle  1,000,000 
nodes.  Thus,  a  volume  of  100  by  100  by  100  nodes  could  be  modeled.  At  1.5  MHz,  this  would 
be  a  5  cm  cube.  It  should  be  noted  that  larger  volumes  can  be  modeled  by  slicing  the  volume  into 
1,000,000  node  rectangles  and  modeling  each  sequentially. 

Within  the  modeled  volume  the  tissue  is  characterized  by  specifying  the  density,  sound 
velocity,  and  acoustic  absorption  at  each  node.  As  a  result  the  interface  is  very  general.  There  is 
much  more  flexibility  in  this  interface  than  is  currently  required  by  the  inverse  optimization 
scheme.  If  we  later  desire  to  optimize  the  sound  velocity  field  within  the  tissue,  or  the  density 
field,  no  modifications  to  the  model  will  be  required. 

Transducer  Model 

If  the  transducer  were  modeled  as  part 
of  the  above  described  volume,  it  would  be 
difficult  to  fit  the  desired  tissue  volume  into 
the  limited  model  volume.  This  problem  is 
addressed  by  solving  for  the  complex  pressure 
field  created  by  the  transducer  being  modeled 
and  applying  this  field  at  one  surface  of  the 
model  volume. 

The  3D  complex  pressure  field  is  first 
found  for  the  transducer  in  water.  This  is 
done  by  matching  a  numerical  solution  to  the 
results  of  calibration  experiments  where  the  temperature  rise  of  epoxy  coated  thermocouples  is 
measured  in  a  water  bath.  Using  this  data,  as  long  as  the  edge  of  the  modeled  volume  extends 
beyond  the  skin  into  the  water  bath,  the  effects  of  the  transducer  can  be  modeled  accurately  by 
finding  the  intersection  between  the  planar  surface  of  the  model  and  the  3D  pressure  field.  Note 
that  a  single  calibration  allows  the  transducer  to  be  positioned  anywhere  with  respect  to  the  model 
volume. 

With  this  user  interface  in  place,  the  model  can  be  treated  as  a  black  box,  where  the 
optimization  scheme  is  ignorant  of  the  details  of  the  model  within  the  box,  but  simply  specifies 
inputs  and  collects  the  output  from  the  model.  Once  the  transducer  has  been  calibrated  (or 
several  transducers),  the  user  simply  specifies  the  transducer  location  and  orientation,  and  the 
nodal  acoustic  properties  that  were  identified  in  Subsystem  1 .  The  model  then  outputs  the 
acoustic  power  deposited  at  each  node  within  the  model  volume. 

Subsystem  4:  Inverse  Technique 

•  Basic  algorithm  development  is  complete.  Optimization  runs  have  been  completed  on 
simple  geometries. 


Complex 
pressure 
specified  on 
this  surface. 
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Inverse  Method  for  Optimizing 
Tissue  Attenuation  Values 

In  the  original  proposal,  it 
was  suggested  that  three  parameters 
could  be  allowed  to  vary  between 
tissue  types.  The  numerical  model 
developed  has  the  capability  to 
change  all  three  parameters. 

However,  in  initial  trials,  it  is 
assumed  that  speed  of  sound,  (c)  and 
impedance,  (Z)  are  known  constants 
for  the  tissues  with  which  we  will  be 
dealing.  Consequently,  attenuation 
is  the  only  parameter  that  changes 
between  tissue  types  and,  thus,  will 
be  the  only  parameter  optimized. 

The  method  discussed  in  the 
proposal  uses  the  Jacobean  matrix  to 
adjust  the  initial  parameter  values  at 
model  nodes  and  iterates  until  the  error  function  is  minimized.  (Here  error  is  some  weighted  sum 
of  the  differences  between  model  predictions  and  experimental  measurements.)  This  method  has 
been  used  in  initial  trials  but  is  being  compared  to  other  methods  to  determine  the  best  method  for 
our  model.  A  second  method,  the  Conjugate  Gradient  method,  is  also  being  evaluated  currently. 

Initial  trials  have  begun  using  each  of  these  methods  to  search  for  minimum  error  points  in 
the  error  function  for  simple  function  spaces.  Figure  16  shows  an  example  where  only  two 
parameters  are  varied.  Given  an  initial  value  for  the  two  parameters,  each  method  iteratively 
searches  for  the  optimum  (minimum)  point  in  the  function.  Essentially  the  only  difference 


Fig  1.  surface  plot  of  function  space  (-.05,  -.03) 
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Figure  16:  This  is  a  surface  plot  of  an  example  error  function 
space.  The  optimization  scheme  is  searching  for  the 
minimum  in  this  function. 


Figure  17:  This  shows  the  rate  of  convergence  for  several  versions  of  the  Conjugate  Gradient  (left)  and  the 
Jacobian  (right)  based  optimization  schemes. 
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between  the  various  optimization  schemes  is  the  efficiency  with  which  they  “slide  down  the  hill” 
to  find  the  minimum  value.  This  efficiency  is  based  on  the  number  of  function  calls  (forward 
model  iterations)  required  to  produce  one  optimization  iteration,  and  on  the  number  of 
optimization  iterations  required  to  reach  the  minimum  error.  Figure  17  shows  trials  using  both 
the  Jacobian  and  the  Conjugate  Gradient  methods  for  optimization. 

Details  of  these  techniques  will  be  discussed  in  the  final  report  and  in  future  papers. 

Concluding  Remarks 

This  report  represents  the  current  state  of  this  project  as  of  September  30, 1997. 
Significant  progress  has  been  made  between  that  time  and  the  delayed  filing  of  this  report.  This 
progress  will  be  discussed  in  the  final  report  and  in  future  papers.  Except  for  the  delayed  filing  of 
this  report,  the  research  progress  is  meeting  the  milestone  schedule  and  we  anticipate  a 
satisfactory  conclusion  of  this  work. 

In  the  4th  year  we  anticipate  the  completion  of  ongoing  work  discussed  above  as  well  as  an 
integration  of  the  efforts  in  the  four  subsections  of  the  project.  Much  of  this  work  will  be 
accomplished  using  agar  phantoms  and  purchased  ultrasound  breast  phantoms.  As  data  becomes 
available  from  live  patients,  it  will  be  used  also.  (Human  patient  or  live  animal  data  will  not  be 
collected  as  part  of  this  research  project,  but  we  will  have  access  to  data  collected  as  part  of  a 
separate  project  supervised  by  Dr.  Roemer.)  As  this  project  concludes,  a  concerted  effort  will  be 
made  to  implement  the  clinical  system  developed  on  the  hyperthermia  treatment  system  under 
development  here  at  the  University  of  Utah.  We  anticipate  the  publication  of  4  or  5  papers 
resulting  from  this  research. 
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